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Discretizing an analytic function on a uniform real-space grid is often done via a straightforward 
collocation method. This is ubiquitous in all areas of computational physics and quantum chem¬ 
istry. An example in Density Functional Theory (DFT) is given by the external potential or the 
pseudo-potential describing the interaction between ions and electrons. The accuracy of the col¬ 
location method used is therefore very important for the reliability of subsequent treatments like 
self-consistent field solutions of the electronic structure problems. By construction, the collocation 
method introduces numerical artifacts typical of real-space treatments, like the so-called egg-box 
error, that may spoil the numerical stability of the description when the real-space grid is too coarse. 

As the external potential is an input of the problem, even a highly precise computational treatment 
cannot cope this inconvenience. We present in this paper a new quadrature scheme that is able to 
exactly preserve the moments of a given analytic function even for large grid spacings, while recon¬ 
ciling with the traditional collocation method when the grid spacing is small enough. In the context 
of real-space electronic structure calculations, we show that this method improves considerably the 
stability of the results for large grid spacings, opening the path towards reliable low-accuracy DFT 
calculations with reduced number of degrees of freedom. 


I. INTRODUCTION 

Real-space grid based techniques are very important 
in disciplines like Quantum Chemistry, Computational 
Physics and Applied Mathematics. A real space ap¬ 
proach is mandatory in the solution of complex problem 
of Partial Differential Equations, as well as for the treat¬ 
ment of complex environments and non-trivial boundary 
conditions. In this framework, the collocation method 
is a straightforward procedure that is used to discretize 
a known function, to express its values in the real-space 
domain. This method represents the most straightfor¬ 
ward and intuitive way to provide numerical coefficients 
to discretize a computational problem. In this method, 
a function / is represented via a set of point values 

fk = f{x k ) , (1) 

where Xk are the sampling points of the simulation do¬ 
main. 

As any discretization method, function collocation in¬ 
troduces an error on the numerical results. This error of 
course decreases while increasing the number of points 
used for the discretization, but its convergence ratio de¬ 
pends of many factors. As an example, imagine that 
the function / represents the charge density of a Poisson 
Equation 

V 2 </> = -4t rf , (2) 

discretized on the points Xk by the collocation method. 
If the multipoles of the coefficients fk do not correspond 
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with the ones of the original function /, a numerical so¬ 
lution of the above equation will not produce a correct 
discretization of the potential <j>, even if the adopted Pois¬ 
son Solver is very accurate. This happens in a number 
of numerical treatment: the discrete multipoles (i.e. the 
momenta of the coefficients fk) of the collocated func¬ 
tions determine the accuracy of the final results. 

In this paper we present a numerical quadrature for¬ 
mula to obtain a set of coefficients fk which can be used 
as a generalization of the collocation method for analytic 
functions. This quadrature scheme is tailored to preserve 
the values of the discrete multipoles of the coefficients, 
to be fixed to the momenta of the original function, even 
for discretization done on grids of large spacings. Such a 
quadrature scheme involves the usage of the Interpolat¬ 
ing Scaling Function (ISF) basis set. 

Similar needs for quadrature formula have already 
been presented in Ref. [I], in the context of the 
grid-point discretization of functions expressed in the 
Daubechies wavelets basis, within the so-called “Magic 
Filter” method. Here we extend and generalize such con¬ 
cept to the real-space discretization of any arbitrary func¬ 
tions. 

In the next section we will illustrate how this problem 
appears in a Real-Space based DFT code, and we show 
the importance of preservation of the monopole in this 
context. We then quantify the discretization errors com¬ 
ing from the collocation method, by explaining our need 
for an alternative formula, and the properties that the 
generalized collocation scheme should satisfy. Then we 
will show the improvements related on the usage of this 
new scheme with supporting results from the BigDFT 
code, showing how the behavior of the results is stabi¬ 
lized for large grid spacings, enabling us to perform low- 
accuracy calculations with reduced number of degrees 
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of freedom. Then we also explain how this quadrature 
method can be used to perform accurately and efficiently 
scalar product between known functions and compactly 
supported basis sets, showing its connections with the 
Magic Filter method. 

II. COLLOCATION AND REAL SPACE 
ELECTRONIC STRUCTURE CODES 

A more elaborated example than the Poisson equation 
presented above is given by electronic structure calcula¬ 
tions on real-space based codes. Density functional the¬ 
ory (DFT) is a widely used and recognized method to 
investigate material properties at the atomistic scale. In 
such ab initio calculations the electronic structure prob¬ 
lem is solved by minimizing the total energy of the system 
which is a functional of the electronic density 

E \p]= T *[p] + \ j p4>\p]dr + E xc [p\ + J Kxt(r)p(r)dr, 

(3) 

where the four terms on the right side of Eq. © are, 
respectively, the kinetic energy, the electrostatic and the 
exchange-correlation energy, and the interaction energy 
with an external potential. Such external potential de¬ 
scribes the interaction of the electrons with the ions, and 
can be considered as the input of the computational prob¬ 
lem. In the Kohn-Sham formalism of DFT, the potential 
Vext appears in the Hamiltonian operator 

'H[ P \ = -\s7 2 + Vks{p\+V^ t , (4) 

where the Kohn-Sham potential Vics[p]( r ) = <j>[p]( r) + 
depends self-consistently on the charge density 
and it is determined during the optimization procedure. 
Its discretization is therefore the primary responsible of 
the quality of the physical and chemical information that 
can be extracted from the DFT calculation. As an ex¬ 
ample, let us now consider a real-space based pseudopo¬ 
tential code. The external potential operator usually is 

bext = Mocal ~b ^nonlocal , (5) 

therefore separated by local and non-local terms. The 
local term, responsible for the electrostatic interaction 
between the ions and the valence electrons, satisfies the 
Poisson Equation 

V 2 Uiocai(r) = -47rp ion (r) , (6) 

where / Oi on (i") is the charge density associated to the 
atomic ions, screened by the core electrons. To have a 
reliable DFT calculation, it is therefore essential that the 
potential V\ oca \ is associated to a p lon with a correct value 
of the monopole, namely the sum of the ionic valence 
charges. 

To illustrate the importance of this consideration, we 
have plotted in Fig. [I] the absolute DFT energy of an or¬ 
ganic molecule, as calculated by the BigDFT code [3] 



FIG. 1. Percentage of the total energy error and the total ionic 
charge of the cinchonidine molecule, C 19 H 22 N 2 O, in function 
of the grid spacing for the collocation method and our pro¬ 
posed one. The total ionic charge of our method has no error 
by construction, which considerably stabilizes the results. 


where V[ oca i is given by norm-conserving GTH-HGH 
pseudopotentials [3], possibly enhanced by a nonlinear 
core correction. The atomic geometry of the molecule 
was optimized using the LDA functional. 

We show that the results become numerically unsta¬ 
ble above the value of 0.7 AU of the grid spacing. It is 
obvious from the figure that this effect is caused, at first 
approximation, by the fact that the total ionic charge is 
not preserved when the grid spacing increases. The dis¬ 
cretization of pi on based on the collocation method limits 
considerably the use of real-space basis sets especially for 
low-accuracy or high-throughput calculations. 

These results are extracted with a particular DFT 
code, but the message is general: each time that the 
collocation method is used, like for example the evalu¬ 
ation of the charge density on a real space mesh (that 
is a fundamental operation, in all DFT codes, needed 
to evaluate E xc [p] and its derivatives), attention must be 
paid to the accuracy of the collocation. Another example 
is the definition of the compensating charge in methods 
like the Projector Augmented Waves ma¬ 
in the following we will present our solution to this 
crucial problem and the potential outcomes for real-space 
based DFT codes. 


III. THE COLLOCATION METHOD AND THE 
INTERPOLATION 

For discretization on uniform grid spacings, the collo¬ 
cation method is well-justified when the original function 
can be reasonably approximated by an interpolation of 
its values on the grid mesh points. Let us consider a one 
dimensional function /. Suppose we want to discretize 
this function on a uniform grid of spacing h and coordi¬ 
nates Xk = hk. Given a family of interpolating functions 
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{Lfe(x)}, if the approximation 

/( x) ~ f L (x) = Y L k (x)f(x k ) (7) 

k 

is reasonably accurate, the collocation method can be ap¬ 
plied. This fact stems from the interpolating property of 
the family {L k (x)}. Indeed, an interpolating family is 
constituted by a set of functions L k , each one associated 
to a point k of the grid, such that L k (j) = S k j- From 
Eq- 0, fL{x k ) = /( x k ) and the continuous representa¬ 
tion of f(x) may be given by /l{x). 

Given the interpolating property, it is also said that an 
interpolating function family is dual to the Dirac deltas. 
In other terms, denoting the above function by the bra- 
ket notations, we have 

I/) - \/l) = Y\ Lk ^ Sk \f) ’ ( 8 ) 

k 

where |dfc) represents the Dirac distribution centered at 
point x k , i.e. (S k \f) = f{x k ). The above defined in¬ 
terpolating property implies that the duality relation 
{6e\L k ) = 6 k £ holds. 

A. Polynomial exactness and discrete multipoles 

The accuracy of the approximation 0 is of great im¬ 
portance for a reliable computational treatment. Clearly, 
such accuracy is intimately related to the family of in¬ 
terpolation functions chosen. The interpolating function 
families are normally constructed using a family of poly¬ 
nomial functions. An interpolating family {L k (x)} is said 
to be of order m if any monomial function x p (denoted 
with |p) in the following), with 0 < p < m is exactly ex¬ 
pressed by the interpolating family, for all x lying within 
a given interval [a, b\. In other terms 

x^Lj(x) = x p ,Vx £ [a, b], 0 < p < m (9) 

j=n a ,n b 

This is the concept of polynomial exactness. Note that 
the index j runs over a set of grid points Xj which might 
lie outside the support [a, 6]. We indicate by [n a ,nf,] the 
minimum interval of grid points needed to obtain the m- 
polynomial exactness in the interval [a, b\. 

The collocation method is therefore meaningful for the 
functions for which the projector operator Y^ k \L k )(5 k \ 
approaches the identity. The polynomial exactness is im¬ 
portant in determining the accuracy of the interpolation: 
a smooth function can reasonably be approximated by its 
Taylor polynomial around a given point. The higher the 
order of the polynomial exactness of the functions L k , 
the better the Taylor expansion of the original function 
would be approximated by the function /^(x), therefore 
the norm of the difference \f — fi\ will be reduced in the 
support [a, 6] by 0(h m ). 

It is easy to understand that this condition is meaning¬ 
ful only when the grid spacing size h is smaller than the 




FIG. 2. Top panel: collocation of a Gaussian function of dif¬ 
ferent standard deviation a on a grid of spacing h = 1. All 
Gaussians are centered in points which lie between two grid 
spacings. It is easy to see that the collocated values are not 
anymore reliable when the ratio h/a grows above 1. This fact 
can also be confirmed by the discrete multipoles of the collo¬ 
cated function, as shown in the bottom panel. It is interesting 
to see that all moments, including the ones which should be 
zero by symmetry, exhibit the same (poor) convergence rate. 

typical oscillations of the function |/) we want to repre¬ 
sent. Evidence is shown in Fig. [2] for a Gaussian func¬ 
tion of standard deviation er, centered between two grid 
points. When a/h < 1, the collocation coefficients f(x k ) 
give a bad representation of the function, as their ampli¬ 
tudes is suppressed everywhere by decreasing a. No in¬ 
terpolating family, even of very high order m, will there¬ 
fore be able to faithfully represent the original function. 
This situation seems unavoidable: as the expansion co¬ 
efficients of the function \f k ) are given in terms of the 
scalar products (Sj \ f ). the grid has to provide a reason¬ 
able sampling of the function / such as to exhibit the 
0{h m ) convergence rate. 

When the grid is small enough, one might think that 
also the discrete multipoles, i.e. 

m p[Il\ = Y xP jf( x ^ ’ 
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( 10 ) 
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follows the same convergence ratio, while approaching 
their exact values given in terms of /. However, it is 
easy to see that this is not the case. The reason is that 
the dual function is always represented by a Dirac distri¬ 
bution, and it is therefore independent of the quality of 
the interpolating family. 

As shown in Figj2j the collocation method gives inac¬ 
curate results for the discrete multipoles when a < h. In 
addition, their convergence ratio is of 0(h), as the coef¬ 
ficients /( x k ) cannot depend on the interpolating func¬ 
tions. At variance from the evaluation of other quantities 
like function derivatives, increasing the order of the inter¬ 
polating function will not change the convergence ratio 
of the discrete multipoles of lower order. 


B. The need of an alternative formula for 
collocation 

Being |/) an analytic function, the accuracy of the mul¬ 
tipoles M p [fi] is a quantitative evaluation of the accu¬ 
racy of /i which is more severe from the one provided 
by the function difference, as it cannot be modified by 
varying the family of interpolating functions. We would 
like to have a different quadrature formula, such that the 
number of preserved discrete multipoles of the original 
function is of the same order of the family of interpolat¬ 
ing functions chosen, regardless of the value of the grid 
spacing. 

As pointed out in the previous section, we need to de¬ 
fine an alternative set of dual functions (Lj |, such that 
the multipoles of the original functions are preserved up 
to order to. In other terms, we search for a family of dual 
functions such that 


cients {fj}, { gj }, the product function should satisfy 


\fg) = J2\ L j)(Lj\fg) 

3 

= ^2 dxLj(x)L k (x)Le(x)f k gt\Lj) 
j,k,e J 

~ £ !,rL<h ■ ( 13 ) 

3 

This is possible only if the dual functions are such that 
J dxl j (x)L k (x)L i (x) = 5 jk Sjt , (14) 


for all the points j , k, £ where the coefficients of the func¬ 
tion are not zero. Note that in the traditional colloca¬ 
tion, when (Lj \ = (<5j|, Eq. (14) is always satisfied due 
to the interpolating property of Lj, i.e. Lj(x k ) = 5j k . 
That proves that the above property is a feature of the 
collocation method. 

We therefore search for a family of bi-orthogonal func¬ 
tions | Lj), (Lj |, such that all the following properties 
hold: 


Collocation coefficients: For dense grid spacings, the 
action of (Lj | should coincide with the Dirac dis¬ 
tribution; 

Biorthogonality: Even though it is not be necessary 
that (Lj\L k ) = Sjk, for any j, k, at least the “weak” 
duality relation 

f Lj(x) L i(x)fi= fj (15) 

i 

should be verified for the points j where fj ^ 0. 


M p [f L ] = '^2,x P j(Lj\f) = (j>\f),0< p < m . (11) 

j 

It is easy to see that this condition can be obtained by 
imposing the TO-polynonrial exactness of the dual set L k 
for all x lying in the support x[f] of the original function. 
In other terms, if the set of L k is such that 


Closure wrt products: The triple product relation of 
Eq. (fl4|) should hold; 


Polynomial exactness: The functions Lj(x) would 
guarantee in this way the nrultipole-preserving 
property. 


In Sec. [A] of the Appendix, we demonstrate that all these 
properties are met for the quadrature formula 


Y^Lj(x)x P j = x p ,0 <p < m,x G x[f] ( 12 ) 
j 

then the multipole preserving property would be guaran¬ 
teed. 

However, we would also like to generalize the colloca¬ 
tion method, rather than to replace it by a new quadra¬ 
ture formula. Firstly, we want to impose that for arbi¬ 
trarily small grid spacings, Lj\f) —► f(xj). Moreover, a 
notable advantage of the collocation method is its closure 
with respect to function products. In other terms, given 
two collocated functions, |/), | g), with collocation coeffi¬ 


fj = (<A) m) |/) = J dxp ( j n \x)f(x) , (16) 

where is a family of ISF of order to. The fj co¬ 


efficients defined by Eq. (16) may therefore be used at 


the place of the function point values f(xj). With this 
choice, we are guaranteed to preserve the first to multi¬ 
poles of / during the discretization procedure. For grid 
spacings which are small enough, we recover the usual 
behaviour of the collocation method thanks to Eq. (A41. 
Fig. [3] provides evidence of this. 

When the coefficients fj are sensibly different from 
f(xj), the values of the coefficients fj should be rather 
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FIG. 3. Top panel: evaluation of quadrature formula of 
Eq. (161 for the same Gaussian functions of Fig. Inter¬ 
estingly, the results for high values of u/h differ considerably 
from the point values: the discretization coefficients may even 
become negative for sharp Gaussians. On the other hand, the 
discrete multipoles M p (bottom panel) of the coefficients agree 
much better with the expected values, and exhibit the correct 
0 {h v ) convergence ratio. 


considered as quadrature terms. However, we might in¬ 
terpret these coefficients as an optimal generalization of 
the collocation method, suitable for grid spacings that 
are larger than the oscillation of functions we would like 
to discretize. This generalization is optimal in the sense 
that the loosening of the accuracy in the multipoles of 
order p, which is unavoidable for large grid spacings, still 
exhibits the correct convergence ratio 0(h p ). 

It is easy to see that, by using a three-dimensional 
separable ISF basis 

= A m) ( x M m) (y)^( z ) > ( 17 ) 

our method can be generalized straightforwardly to 
three-dimensional grids, especially for separable func¬ 
tions. In the following section we will illustrate the ad¬ 
vantage of this method for electronic structure calcula¬ 
tions. 


IV. CONSEQUENCES IN ELECTRONIC 
STRUCTURE CALCULATIONS 

We illustrate our idea using the BigDFT code [2j which 
is based on Daubechies wavelets to express the electronic 
wavefunctions and on interpolating scaling functions for 
the electronic density and the Kohn-Sham potential. All 
tests in this article will be done with the LDA functional. 
We have checked that these results are the same for differ¬ 
ent functionals and also for the Hartree-Fock approach. 

The BigDFT code has an adaptive mesh with one level 
of refinement and the corresponding parameter hgrid 
specifies the grid spacing of the coarse resolution. The 
finer resolution which is only used near the nucleus so has 
a twice finer grid step by construction. As mentioned in 
Sec.[n| norm-conserving GTH-HGH pseudopotentials 0 
are used in the BigDFT code. They are built with Gaus¬ 
sian functions for p lon and V n oniocai- Using the colloca¬ 
tion method, BigDFT needs a grid step of the order of 
the standard deviation er parameter of the Gaussian func¬ 
tion. As an example, for the case of the hydrogen atom, 
a has a value of 0.2 atomic units that obliges to use a 
grid spacing of the same value i.e. a/h > 1. In the case 
of BigDFT, this means that the input parameter hgrid 
should be of the order of 0.4 AU in order to have the 
finer mesh of the same resolution as the Gaussian stan¬ 
dard deviation parameter. 


A. Accuracy in absolute energy 


In Fig. |4j we show the percentage of the difference 
of the total energy from the reference calculation with 
hgrid= 0.15 in function of the grid spacing. The quadra¬ 
ture formula (16) has been used with ISF of order m = 


16. For a grid step greater than 0.6, the collocation 
method gives an error bigger than 1%, drastically in¬ 
creasing. On contrary, our multipole preserving method 
is more stable given an error of 1% for hgrid=0.9 with 
an accuracy of two orders of magnitude up to a grid step 
of 2 which corresponds to five times the a value of the 
Gaussian function. 

Another intrinsic artifact of the real-space methods, 
especially the finite difference scheme or any methods 
which uses collocation technique, is the egg-box er¬ 
ror mm- The discretization procedure is not invari¬ 
ant by global translations and rigid rotations. So dif¬ 
ferent relative positions of the ions with respect to the 
mesh of the simulation domain change the discretized 
values of V e xt- This egg-box effect can be considerably 
reduced by our method, for large grid spacings. In the 
figure [5j we have plotted the maximum variation of the 
absolute energy of a H 2 molecule when rotating its main 
axis, while keeping fixed interatomic distance. As the 
relative positions of the centers of the pseudopotentials 
vary strongly with the molecule orientation, this is a 
good (even severe) estimator of the egg-box error. We 
can see that the multipole preserving method decreases 
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Grid step (A.U.) 

FIG. 4. Percentage of the error on the H atom total energy 
using the LDA functional in function of the grid step for the 
collocation and the multipole preserving methods. The refer¬ 
ence total energy is calculated with a grid spacing of 0.15. 


Number of basis set functions 
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FIG. 5. Maximal egg-box error in meV per atom for any 
rotation of the Ha molecule in function of the grid step. 


considerably the egg-box error in an intermediate range 
between 0.4 and 0.9 where the egg-box energy becomes 
bigger than 100 meV/atom. Not surprisingly, the tra¬ 
ditional collocation method totally fails to capture the 
rotation of the H2 molecule with a very large error. 

We would like to point out that the quasi-variational 
behaviour, typical of any real-space method based on 
analytic potentials, is greatly improved for large grid 
spacings. This point is interesting because it is possible 
to reduce the degrees of freedom by increasing the grid 
spacings, therefore decreasing the accuracy, but without 
spoiling the physical meaning of the results. This could 
be an interesting alternative to tight-binding methods 
based on DFT which need fine-tuning of many parame¬ 
ters to represent properly a molecule with an accuracy 
of the order of 100 meV per atom. In the upper x axis 
on the same figure [5j we indicate the number of basis set 
functions i.e. the number of degree of freedoms used for 
each grid step. These numbers vary as the inverse of the 
cube of the grid step and the CPU time is directly propor¬ 
tional to the number of the basis set functions. Thus it is 
very advantageous to increase the grid step parameters to 
save a lot of CPU time if the method is numerically stable 



FIG. 6. Bond distance and error of the relaxed H2 molecule 
in function of the grid spacing. 


especially when we explore the potential Energy Surface 
(PES) with some methods as Minima Hopping [14| or the 
Activation Relaxation Technique m- 


B. Accuracy in Geometry optimization 

To illustrate this idea to have a reliable PES, we have 
optimized, in the figure [6] the H 2 molecule in function of 
the grid spacing using the FIRE [16] method because this 
method is robust when the egg-box error becomes large. 
Below a value of 0.8 for the grid spacing which corre¬ 
sponds to 2 times the standard deviation value of the 
Gaussian of the pseudopotential, the collocation method 
still provides manageable results, giving an error less 
than 0.1 angstroem for the bond length. For a larger 
value of the grid spacing, the collocation method becomes 
strongly unstable. On contrary, our method is stable 
numerically on a wider range, even though, clearly, the 
precision of the results decreases. 

For an exploration of the PES with an accuracy of the 
order of 0.1 angstroem a grid spacing of the order of 1 
AU permits to decrease the number of basis set function 
by two order of magnitudes (see Fig. |5j) accelerating in 
the same range. Such low-accuracy setup might be suffi¬ 
cient for pre-screening calculation in workflows like high 
throughput funnels. In the same spirit as high through¬ 
put screening, using the same input parameters and the 
same code, varying only the grid spacing i.e. the accu¬ 
racy, could constitute an interesting approach to have a 
preliminary investigation of interesting minima and sad¬ 
dle points. 

Finally, the multipole-preserving scheme permits to 
use larger grids space to calculate organic molecules con¬ 
taining hydrogen but also carbon, nitrogen and oxygen 
atoms where the coefficients of the pseudopotentials are 
between 0.2 and 0.35. In Fig. [7]and Fig. |I] we show that 
the total energy of the cinchonidine molecule is numeri¬ 
cally stable in function of the grid step. 
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support, ranging from i a to % = i a + n , and its analytic 
moments (p\4>), the coefficients 

n— 1 

( 2 0 ) 

9=0 

might help us in defining a “smoothed” function 

l^ (n) > = X>f Vl n) > ( 21 ) 

i—i a 

such that 


FIG. 7. Percentage of error of the total energy of the cin- 
chonidine C 1 QH 22 N 2 O molecule for the collocation and the 
multipole-preserving schemes. 


V. COLLOCATED SMOOTHENING 
PROCEDURE FOR A COMPACTLY 
SUPPORTED FUNCTION 


(0 (n) l/) = EX*** = W)+0(h n ) , (22) 


where the /* coefficients are calculated by Eq. (16) and 


can be replaced by f(x{) only when h is small enough. By 
construction, Eq. ( [22] ) is exact for polynomial functions 
of order less than n. 


In this admittedly more technical section we show 
that our quadrature formula might be used as a reli¬ 
able approach to calculate the scalar product between 
our original function / and a compactly supported func¬ 
tion (j>. This constitutes an important consequence of our 
method when known functions have to be represented in 
compactly-supported basis sets. 


A. Lagrange polynomial as the direct basis 

Let us now come back to the approximated function 
Jl defined in 0. On a support ranging on a interval 
[ia, ib = ia+n], in Sec. [A] of the Appendix, we show that 
a reliable approximation of / can be given by defined as 


B. The “Magic Filter” method in BigDFT code 


The above result has been already pointed out in the 
context of the so-called “Magic-Filter” method, provid¬ 
ing an optimal set of quadrature coefficients to convert a 
function expressed in the Daubechies wavelets basis | </> M ) 
into its point values on a uniform grid. In Sec. [B] of the 
Appendix we show that the usual procedure for calculat¬ 


ing these filters leads to equation (20). 


Such method is very useful in the BigDFT code to 
express efficiently the local potential matrix elements in 
Daubechies wavelet basis. In the BigDFT code, a Kohn- 
Sham orbital | \B) is expressed in Daubechies wavelets 


w = E c ^i < w c M = (^i^) 


(23) 


h(x) = ^2,C { - a ' Xb \x) J ip ( J m \x)f{x)dx , (18) 


where 




(x) are interpolating Lagrange polynomials, 


defined as in Eq. (A7) 


These properties are interesting any time we have to 
perform the scalar product of the approximated function 
/l with a compactly supported function |</>), of support 
in the interval \i a ,ib]- Indeed, the integral 


Wl) = = f>|Tp ! ’ ) >(<^ ) | f)+0(h n ) , 

(19) 

will be exact if | /) is a polynomial function of order less 
than n. This happens because also the Lagrange poly¬ 
nomials £^ a ' lb ' ) satisfy n-polynomial exactness in [i a ,ib] 
interval. 

The above results, can be also provided in a more gen¬ 
eral sense. Given an arbitrary function | <j>) of compact 


As pointed out before, the Daubechies wavelet basis is 
a compact support orthogonal basis |0 M ) able to exactly 
express the polynomials up to a give order m. However, 
the basis functions 

Mx) = “ h) ( 24 ) 


have the peculiar property of being less smooth than an 
order ?n-polynomial. Therefore, the naive collocated val¬ 
ues ^(x^) would not provide an efficient discretization. 

We can therefore use Eq. (21) to define an improved 
collocation method. We have seen that it is more accu¬ 
rate to express the Kohn-Sham orbitals in real space by 
the following smoother function: 


®l(x) = E c ^^ 2m) ( a; ) = wf 51 55 W m-Z (x) 

n V h 11 i=l—m 

, m 

= - 7 ^ E ^ c i-3 w wf m \ x ) » (25) 

v 2 = 1 —m 7 








1.2 



t 


the Daubechies wavelet basis and a real-space descrip¬ 
tion in a generalized collocation scheme. As shown in 
mm Eqs. ( f26j ) and ( [28]) show that this passage ma¬ 
trix is unitary up to (D(h Im ). If the Kohn-Sham orbital 
would be analytically known, we could apply Eq. (22), 
and our generalized collocation scheme would reduce to 
Magic Filter method. For this reason, even though it 
is generally applicable, the multipole-preserving quadra¬ 
ture presented in this paper perfectly reconciles with a 
Daubechies wavelets computational treatment equipped 
with the Magic Filter method. 


VI. CONCLUSION 


FIG. 8. Smoothed Daubechies scaling function cp^ 16 \t), with 
the prescription given by Eq. (211. The values at the integer 
points Wi = are the “magic filter” coefficient for the 

original Daubechies scaling function of order m — 8, indicated 
by 


whose collocated values f(xj ) would preserve the multi¬ 
poles of that can be derived from the original expres¬ 
sion in terms of coefficients. This function is expressed 
in terms of the smoothed Daubechies scaling functions 
cf>\} 6 \t), plotted in Fig. [sj Therefore, when starting from 
the coefficients c^, the real space values of ip might be 
given by the formula 


fi 


, m 

V ct - 


Wi 


(26) 


The results of this paper might be useful to define the 
inverse relation. Given a set of real-space point values 
fj, these coefficients might be interpreted as generalized 
collocation values. With this interpretation we are able 
to write the piecewise polynomial expansion of the Kohn- 
Sham orbital valid on a interval of size 2m around a grid 
point p,: 


$l(i) 




j—l—m 


(27) 

This piecewise polynomial function would have the ex¬ 
pansion coefficients in Daubechies wavelets basis given 
by the equation 




1 m /*m 

j=l — m 


= Vh ^2 . (28) 

j=l-m 


This result shows that the “Magic Filter” method 
can be seen as the optimal passage matrix between 


Collocation method is a universally applicable pre¬ 
scription for the numerical discretization of functions. 
However it suffers from an intrinsic limitation: highly 
oscillating functions cannot be well represented on a grid 
if the spacing is too large with respect to the typical 
length of the oscillations. Therefore the accuracy of the 
collocation is rapidly spoiled as soon as the grid spac¬ 
ing becomes too large. Unstable results might occur if 
the numerical implementation is done in such a regime. 
This limitation implies that there is a upper limit for the 
grid spacing in a real space based DFT code, and con¬ 
sequently a lower limit for the number of computational 
degrees of freedom. Results become rapidly meaningless 
when these limits are overcome. 

With this paper, we have presented a method to gen¬ 
eralize the collocation of arbitrary analytic function on 
large grid spacing, without spoiling the accuracy of the 
discretization. The collocation values might be replaced 
by the scalar products of the analytic function with the 
basis of Interpolating Scaling Functions. For analytic 
and separable functions like Gaussians, this prescription 
is very simple and easy to implement in three dimensions 
(see e.g. 0), and tends to the point values when the grid 
is fine enough. This method has been implemented in the 
BigDFT code, which uses a real-space based description 
using Daubechies wavelets. Thanks to the inclusion of 
this method, the code exhibits numerical stability over 
a wide range of grid spacings, not accessible with tradi¬ 
tional collocation. However, the implementation of this 
method is unrelated to Daubechies wavelet basis set, and 
can be used in other DFT codes and even in different 
contexts, like for example the definition of compensating 
charges in Fast Multipole Methods. 

Having said that, we would like to point out that 
our method is focused in providing stabilisation of low- 
accuracy results, that would be unaccessible if traditional 
collocation is used. If a high-accuracy calculation is 
needed, the grid spacing has to be adjusted in the re¬ 
quired range. Nonetheless, as it preserves numerical sta¬ 
bility, our multipole-preserving approach is guaranteed 
to be better than traditional collocation. 

The outcomes of this method are very important, as it 
enables us to use larger grid spacings even for hard pseu- 
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dopotentials and to perform coarse-grained DFT calcu¬ 
lations. As the large grid behaviour of the code is highly 
stabilized with this method, the user is now able to per¬ 
form low-accuracy DFT calculations with reduced num¬ 
ber of degrees of freedom. This is fundamental in view 
of rapid exploration of the energetic features of a system 
at DFT level, or to accelerate the convergence of iter¬ 
ative molecular calcualtions m • Notable examples are 
the Potential Energy Surface explorations of systems at 
the nanoscale, as well as the recently established held of 
high-throughput calculations for material design. 
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FIG. 9. Plots of interpolating scaling function and 

corresponding lifted wavelet with 16 vanishing moments. 


Appendix A: Interpolating scaling functions 

Interpolating scaling functions (ISF) 0 arise in the 
framework of wavelet theory 0 El- They are one¬ 
dimensional functions, and their main properties are: 

• The full basis set can be obtained from all the trans¬ 
lations by a certain grid spacing h of the mother 
function <centered at the origin.We indicate 
the basis set with 

p^\ x ) = ^(^-i) (Al) 

• The mother function p( m l is symmetric, with com¬ 
pact support from —m +1 to m — 1. It satisfies the 
interpolating property p^ m \j) = Sj. 

• They satisfy the refinement relation 

m— 1 

v {m \x)= Y W m) (2 X~j) (A2) 

j =- m -\-1 


The ISF families exhibit polynomial exactness: indeed, 
the so-called lifting procedure allows to define a set of 
functions {4’j} the lifted wavelets - that are both or¬ 
thogonal to pf ] and have (at most) m vanishing mo¬ 
ments. As the multi-resolution basis formed by the pj 
at lowest resolution and the lifted wavelets ifj at all the 
resolution levels forms a complete set, this proves that 
the basis of the pj exhibits polynomial exactness up to 
order m. Fig. [9] shows an interpolating scaling function 
pl ie l (t), together with the corresponding lower resolution 
lifted wavelets. 

The polynomial exactness of the pj basis, together 
with their compact support, allows us to demonstrate 
the collocation property: we can demonstrate that 

lirn [ d xpf l \x)f(x) = lirn f df<p (m) (t-j)f(ht) = f(hj) , 
h —>0 J J h —>0 J 

.(A4) 

where we have expressed x = ht in terms of the dimen¬ 
sionless unit t. To show this, let us suppose that the 
function f{x) can be well approximated by its Taylor 
polynomial close to the point x :) = hj: 


where the hj ’s are the elements of a filter that char¬ 
acterizes the wavelet family, equal to p m (j/2) in 
the case of ISF, and m is the order of the scaling 
function. Eq. (A21 establishes a relation between 


the scaling functions on a grid with grid spacing 
h and the ones defined higher resolution level with 
spacing h/2. 


The filters in Eq. (A2) are defined such that (as 
proven in Ref. 0 ) that the lowest m moments of 
the scaling function are all vanishing but the first, 


i.e.: 


(*><”‘V>=/^Kd* = V 0 < p < m 


(A3) 


This enables us to show that (see e.g. 0) 
(<p- m) |p) = x p : 


m —1 


m) = Y - jy + o{h m y ,, 


g=0 




which would lead, thanks to Eq. (A31, to 
f d tp^\t-j)f{hf) 


= y f d t{typ {m) {t)+ 

9=0 q ' 


(A5) 


+ 0(h m ) = f(hj) + 0(h m ) , (A6) 


which proves Eq. (A41. 

We have seen that for a grid spacing which is small 
enough, the coefficients fj = (p^lf) approach the col¬ 
location coefficients f(xj). Therefore in this case any set 
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of interpolating functions with polynomial exactness of 
order m would be a good choice for Lj, as the corre¬ 
sponding f L would be similar up to 0(h m ). As a matter 
of fact, ISF basis can be treated as an orthogonal basis, 
even this is not true: non-orthogonality effects are visible 
only up to 0(h m ). 


However, even in this case, it might be useful to iden¬ 
tify the direct basis Lj which better generalizes the col¬ 


location approach. The Lagrange polynomials 

ib(i^k) _ . i b -i a -1 

ct ib \t)= n = E 44 9 ’ ke[i a ,i b ] 

i=i a q =0 

. (A7) 

might constitute a basis of interest. The matrix A l £’ bb is 
the inverse of the Vandermonde matrix, i.e. 

Y A )\q b i q = Sij Vi, j G [i a ,ib\ ■ (A8) 

q =o 

In the following we show by using this basis set as the 
direct interpolating basis, we recover the biorthogonality 
and the triple product property (Eq. (HD- 

Let us consider a set of discretization points ranging 
from i a to i b = i a + n. The properties of ISF imply 


71—1 

(^ ra) 4' 0 '“ ) ) = ^144* = 5lk G [ i a,ia+n],n < m , 

9=0 

, 71—1 71—1 

dx Xl4 T x< i = Si ^ ik e [ia,ia + n],n < m/2 . 

p —0 q —0 


(A9) 

(A10) 


Therefore, for all the i lying in the interval [i a , i b = i a +n], 
the basis set constitutes a biorthogonal 

basis generalizing the collocation method to 0{h n )[? ]. 


Appendix B: Magic Filter, the original idea 


Given a set of 2m-family Daubechies scaling functions 
c t> centered on a uniform mesh of spacing h , the expansion 
coefficients of a given function f(x) in this set are defined 
as 


= 4 y - n)f(x) = Vh j d t<j>(x - n)f(ht) . 

(Bl) 

This discretization of the function / in this basis set has 
an algebraic h m convergence rate. In other terms 

f(x) - 4 £ c ^(| — (a) = 0(h m ) . (B2) 


A wavelet quadrature in this context is based on the 
idea of approximating the expression above by a collo¬ 
cation formula. In other words, we should define some 


coefficients ru,, where i = 1 — m, ■ • • ,m such that 

= VhY t w j - li f(hj) + VhO{h 2m ) . (B3) 

i 

This can be possible only if the above formula would 
give the ex act r esult whe n f(x)=x p ,0<p< 2to. By 
comparing (Bl) and (B3) in the case of polynomials we 
thus find the equation defining the Magic Filter: 

y Wj ( j + i) p = / dxcj)(x)(x + i) p V*, 0 < p < 2?n , 

( B4) 

which is solved V* if and only if 

Wjj p = J d x(j){x)x p = M p VO < p < 2to , (B5) 


the Magic Filters are given by Eq.(10) of Neelov and 
Goedecker paper [T] : 

m r. 

iv k = = J dxc/)(x)C { k ~ m,m) ( x ) ■ 

7 = 1 — 771 

. ( B6 ) 

This equation shows that the magic filters can be viewed 
as the expansion coefficients of the Lagrange polynomials 
in the Daubechies scaling functions basis, and is identical 
to Eq. (20). Also the paper from Johnson [5] can be used 
as a reference in this regard. 
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